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SUMMARY 


For  a  stationary  flow  in  a  cylindrical  container  of  the 
Couette  type  in  an  outer  radial  zone,  and  of  zero  velocity  in  an  inner 
radial  zone,  the  normal  mode  equations  are  derived.  For  negative  wave 
numbers  in  the  8-direction,  these  equations  are  found  to  have  a  simi¬ 
larity. 


The  eigenvalues  are  calculated  by  initial  value  methods  em¬ 
ploying  the  Runge-Kutta-Gill  integration  procedure.  Values  of  the 
dependent  function  and  their  derivatives  at  the  singularity  are  calcu¬ 
lated  by  linear  extrapolation  coupled  with  continuity  requirements. 

Tables  of  eigenvalues  for  various  slenderness  ratios  of  the 
cylinder  and  various  radial  nodes  are  given  for  0-wave  numbers  of  -1. 
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INTRODUCTION 

The  object  of  the  following  investigation  is  to  determine  the 
eigenfreauencies  of  a  viscous  fluid,  contained  in  a  rotating  cylinder, 
during  the  initial  spin-up  period.  In  order  to  make  the  problem  trac¬ 
table,  an  inviscid  fluid  is  assumed  and  the  fluid  is  divided  into  two 
zones;  an  outer  radial  zone  in  which  a  Couette  type  flow  prevails  and 
an  interior  zone  which  is  static.  The  initial  flow  distribution,  where 
the  interface  between  the  two  zones  is  assumed  to  decrease  from  the 
outer  radius  to  zero  in  a  quasi-static  manner,  is  assumed  to  be  time  in¬ 
dependent  . 


FORMULATION  OF  THE  EIGENVALUE  PROBLEM 

In  the  following  we  are  interested  in  the  characteristic  fre¬ 
quencies  of  an  inviscid  liquid  contained  in  a  cylindrical  container. 

The  liquid  is  assumed  to  have  an  initial  stable  motion  of  the  Couette 
type  for  an  outer  radial  zone  whereas  the  interior  is  assumed  to  be  at 
zero  velocity.  Thus,  let  the  container  be  of  radius  a  and  height  2c. 
Then  the  stationary  flow  in  the  outer  zone  is  given  by 

u«0  a>r>b 

o  —  — 

vq  »  au>[(r/a)  -  (e^a/r)]/(l  -  e^)  0  _<[e  *  b/a]  _<  1 

w  *  0 
o 

In  the  inner  zone, 
ui  «  0 

v1«0  b  >  r  >  0 

w^  ■  0 

where  b  is  the  radius  of  the  interface  between  the  static  interior  and 
the  moving  exterior;  u,  v,  and  w  are  the  radial,  longitudinal  and  axial 
components  of  velocity  in  a  cylindrical  coordinate  system  with  the 
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z-axis  aligned  along  the  axis  of  the  cylinder.  The  origin  of  the  coor¬ 
dinate  system  is  located  at  the  bottom  of  the  cylinder  so  that  the  z- 
coordinates  of  the  end  faces  are  given  by  0  and  2c.  The  velocity 
components  in  the  inner  and  outer  zones  are  distinguished  by  the  sub¬ 
scripts  i  and  o  respectively. 

For  convenience,  we  consider  the  two  zones,  namely,  the  exterior 
zone  of  Couette  flow,  and  the  interior  zone  of  zero  velocity,  separately. 
These  are  then  coupled  by  boundary  conditions  imposed  at  the  interface, 
where  the  radial  velocities  and  pressures  corresponding  to  the  two  regions 
are  required  to  be  equal.  Considering  the  exterior  zone  first,  the  con¬ 
tinuity  and  momentum  equations  with  the  appropriate  boundary  conditions 
are: 


3u 

_ c 

at 

3v 


+  (o  -V)  u  -  vz/r  =  -  f  -  &) 
K  o  J  o  o  3r  '‘D  ' 


t— ^  +  (o  *v)  v  +  u  v  /r  =  -  —  (— ) 
3t  ^  o  ;  o  oo  r  30  } 


3w 

_ c 

3t 


+  (o*v)w  =  -  (— ) 

v  o  ;  o  3z  '‘p  J 


Q.y  =  u  —  +  SL  il_  4.  w  JL_ 

■o  Uo  3r  r  36  o  3z 

u  3u  ,  3v  3w 

+_0  +I_£ +  _o  =  0 

r  3r  r  36  3t 


The  boundary  conditions  are 


r  -  a;  u  »  0 
0 

z  ■  0,  2c;  w  =  0 
o 


Moreover,  at  the  interface,  we  have 


“of  "  Uif;  Pof  "  Pif 


where  subscript  f  refers  to  the  interface. 


(3) 


(4) 
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where  is  the  equilibrium  pressure  at  the  interface. 

We  now  assume  perturbations  due  to  a  radial  displacement,  n>  of  the  in¬ 
terface  from  its  equilibrium  position  corresponding  to  r  =  e. 


u 

e 

u' 

0 

0 

V 

B 

V  +  v* 

0 

0 

w 

B 

w* 

0 

0 

p 

frvi£ 

0 

a) 

dy  +  P'  +  P- 
o  1 

Substituting  the  perturbations  into  the  dimensionless  Euler* s  equations 
of  motion,  and  neglecting  second  and  higher  order  quantities, 


(ID 


3u*  „  3u* 

o  ,  V  _ o 

3t  r  36 


2V 


3P* 

_ o 

3r 


3v* 

It 


o  ,  dV  V  3vo 

+  U  —  +  — 


u*V 

o 


o  dr  r  36 


-  3P  * 
1_ _ o 

r  36 


3w*  ..  3w*  3P* 

_ o  V  o  _  o 

3t  r  36  3z 

u*  3u *  .  3v*  3w* 

_o  +  _q  +  I_o+_£ 

r  3r  r  36  3z 


The  boundary  conditions  become 


(12) 


r  -  1;  u’  -  0 
o 

z  -  0,  2c/a;  =  0  (13) 

The  boundary  conditions  due  to  displacement  of  the  interface  are  obtained 
in  the  following  manner: 

Let  the  interface  be  at 


r  -  e  +  n  (14) 

Substituting  into  the  last  of  Equations  11  and  neglecting  second  and  higher 
order  terms, 


pi  +  p;£ 


05) 
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At  the  interface,  we  have  from  14 


u 


la 

o  3t 

where  second  and  higher  order  quantities  in  the  perturbations  have  been 


(16) 


ignored. 

We  now  consider  the  interior  zone.  The  Euler's  equations  of 
motion  and  the  equation  of  continuity  are  once  again  given  by  Equation  3 
in  a  non-dimensional  form.  The  boundary  conditions  are 


z  ■  0,  2c/a;  w. 


0 


(17) 


r  ■  0;  u^,  v^  and  w^  are  bounded. 


At  the  interface,  they  are 


u 


if 


u 


of 


rif  of 

We  once  again  assume  perturbations  due  to  radial  displacement,  q,  of  the 
interface  from  its  equilibrium  position  corresponding  to  r  8  e. 


(18) 


u. 


w, 


W 


p  «  p*  +  p 
i  i  rl 


Substituting  in  Equation  3, 

3u^ 


3t 

Hi 

3t 

!! 
<*t 
u’ 


3P^ 

17" 
i  . 

r  36 

3z 


1  3Pi 


(19) 


_i  S  1^1  K 

r  +  3r  +  r  30  +  3z 


(20) 
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The  boundary  conditions  become 


z 

r 

At  the  interface 


0,  2c/a,  w^  *  0 

0;  u^,  v|  and  are  bounded. 


la 

at 


NORMAL  MODES 


(21) 


(22) 


In  accordance  with  the  usual  procedure  of  treating  characteristic  value 
problems,  the  perturbations  are  analyzed  into  normal  modes.  In  view  of 
the  boundary  conditions,  it  is  natural  to  suppose  that  the  perturbations 
are  given  by  quantities  which  have  a  (r,  0,  z,  t)  dependence  given  by 

i(K  t  +  m  0) 

u*  ■  U  (r)  cos  [h  7raz/2c]e  u 
oo  o 

i(K  t  +  m  6) 

v*  ■  B  (r)  cos  [h  TTaz/2c]e 
oo  o 

i(K  t  +  m  0) 

w*  *  W  (r)  sin  [h  7raz/2c]e  1 

o  o  o 

i(K  t  +  m  0) 

P*  *  R  (r)  cos  [h  Traz/2c]e  0  0  (23) 

o  o  o 

Where  Kq  is  a  constant  (which  can  be  complex) ,  m^  is  an  integer  (which 

can  be  positive,  zero,  or  negative)  and  hQ  is  the  wave  number  in  the  z- 

direction  (which  can  be  1,  2,  3,  etc.).  U  ,  B  ,  W  and  R  are  functions 

ooo  o 

or  r  only. 

In  an  analogous  manner,  for  the  interior  zone,  the  perturba¬ 
tion  quantities  may  be  assumed  to  be  of  the  form 
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i(K  t  +  m  6) 

u*  *  U.  (r)  cos  (h  iraz/2c)e 

i(K.t  +  m  6) 

v'  ■  B  (r)  cos  (h  iraz/2c)e 

i(K.  t  +  m.e) 

w*  *  W,(r)  sin  (h  7raz/2c)e 

i(K,t  +  m.e) 

*  R^(r)  cos  (hiiraz/2c)e  i 

where  K^,  m^,  and  are  defined  as  before.  U^,  B^,  W^,  and  are  once 
again  assumed  to  be  functions  of  the  radius  r  only. 

Substituting  Equation  23  into  Equation  12,  we  obtain 


i(K  +  m  -)u  -  —  B 
v  o  o  r'  o  r  o 


(^  +  -)u  +  i  (k  +  m  -)b  -  -i  m  — 

vdr  r'o  ^  o  o  r'  o  or 

„  h  ira 

(k  +  m  -)w  -  -  -2—  i  R 
v  o  o  rJ  o  2c  o 


dU  U  \  B  h  na 

“ —  +  —  +  i  m  —  +  — " —  W  =0 

dr  r  /  or  2c  o 


obtain 


In  a  similar  manner,  substituting  Equation  24  into  20,  we 


i  dR± 

U  *  —  - — - 
i  K±  dr 

1  mi 

Bi  =  -7K^Ri 

i  hi1ra 

Ri 

U  dU  im  h  ira 

—  +  ^  +  —  B4  +  ~ —  W  =  0 

r  dr  r  i  2c  i 


(26) 


Since 


u’  =  u' 
of  if 
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we  have 

i(K  t  +  m0) 

U  (e  +  n)  cos  (h.Traz/2c)e 

i(K  t  +  M.) 

-  U  (e  +  n)  cos  (h  uaz/2c)e  0  0  I 

o  o 

Inasmuch  as  0,  z  and  t  are  Independent  variables,  we  require  that  K^, 
and  h^  be  equal  to  Kq,  mQ  and  hQ  respectively.  Thus,  the  subscripts  for 
K,  m  and  h  will  be  omitted  in  what  follows.  Equation  25  can  be  simpli¬ 
fied  to  the  following 

[dU  U  1  0  I  2  2,  2  2™1 

o,  o  m  TT  d  m  ,  ir  h  a  .  ~ 

3r+7-J-7".3?(“)-  b  +  — r  K 


d2U  -  &  U  f-  («r2) 
o  r  o  dr  1  ; 


R  dR 

217m  — “  +  a  ~ — “  i 
r  dr 


where 


n  =  V/r 


a  =  K  +  mft 


In  a  similar  manner,  Equation  26  becomes 


U.  dU4  .  2  2,2  2 
_i.  _ t  _  JL  B  .  h  a 

r  +  dr  "  K  2  +  .2  Ri 

r  4c 

Ui  K  dr  ' 

The  boundary  conditions  are 

UQ(1)  -  0  < 

U^(o)  and  R^(0)  are  bounded 

The  interface  conditions  are  obtained  as  follows:  Since  P^  is  equal  to 
P^,  from  Equations  15  and  22, 
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Substituting  for  and  Pj^  from  Equations  23  and  24  respectively, 

\  /U  y0  *  i(Kt  +  m0) 

R^e)  cos  (hTTaz/2cye 


R  (e)  cos  (h7raz/2c)e 
o 


i(Kt  +  m0) 


(34) 


where  we  consider  n  <<  e. 


fu.  dU  |  fdU  U  |  ,  0 

K{ri  +  ar}-0{ar  +  r}-7uofc^2)Jre 

UQ(e)  -  Ui(e) 

Thus,  the  characteristic  value  problem  consists  of  Equations  28  and  31 
with  boundary  conditions  32,  35,  and  36.  Since  we  are  interested  only 
in  those  motions  which  produce  a  couple  on  the  casing,  interest  centers 
around  values  of  h  =  (2j  +  1),  where  j  is  an  integer. 


(35) 

(36) 


Axisymmetric  Cases 

Setting  m  =  0  reduces  Equations  28  and  31  to  the  following 


d2u 


dr 

d2U. 


i  dU 

0,1  _ o 

2  rdr 


1  .  „2  N  2ft  d  rn  2> 

“  +  N  "7~d7  ) 

r  K 


U  =  0 
o 


dr 


2  2  2  2  2 
where  N  =  tt  h  a  /4c 

The  boundary  conditions  become 


U0(l)  =  0 

U  (e)  =  U.(e); 
o  i 

1^(0)  -  0 


dUo(e)  dUi(e) 


dr 


dr 


(37) 

(38) 


(39) 


* 

The  axisymmetric  case  is  of  academic  interest  only. 


-  9  - 


IHTth 


THE  FRANKLIN  INSTITUTE  RESEARCH  LABORATORIES 


F-B2294 


Case  1,  The  Narrow  Gap  Approximation  for  m  =  0 

In  specifying  ft,  interest  is  in  those  velocity  distributions 
which  are  realizable  in  a  viscous  fluid.  Equations  9  and  29  determine 
the  case  of  primary  interest  as 

0  -  (l  -  e2/r2)/(l  -  e2)  (AO) 

However,  as  e  1,  an  important  simplification  in  the  characteristic 
value  problem  given  by  Equations  37,  39  and  40  is  possible  provided  that 
(1  -  e)  is  small  when  compared  with  (1  +  e)/2.  In  this  case,  Equation  40 
may  be  expressed  as 

ft  =  s  (41) 

s  «  (r  -  e)/(l  -  e)  (42) 


The  first  of  Equation  37  may  be  expressed  as 

,2 


d  Id  1  **2 

71  *7d7-~ -N 

dr  r 


U  =  -  20r  ^  +  AO2  U 

0  k2  L  dr  J  ° 


Using  the  transformation  defined  by  Equation  42,  consistent  with  the 
"narrow  gap"  approximation  of  small  (1  -  e) ,  equation  43  becomes 

2 


ds 


2  ~  al 


it  1  2e 

U  - - -  •  - -  s  U 

o  2  1-e  o 


where 


2  _  ,  _  v  2.  2 

a^  =  (1  -  e)  N 

By  means  of  the  following  transformation, 

2/3  ,2/3  1 

x-*1  b  b'-J 

where 


(A3) 


(AA) 


(A5) 


(A6) 


b2  =  2e/K2(l  -  e) 


(A7) 
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Equation  44  may  be  expressed  as 


x  U 

o 


0 


(48) 


The  general  solution  of  Equation  48  may  be  given  in  terms  of 
either  Bessel  Functions  of  order  1/3  or  somewhat  more  conveniently  in 
terms  of  Airy's  functions  AI(x)  and  BI(x);  thus, 


U  =  A  AI(x)  +  B  BI(x) 
o  o  o 

where  Aq  and  Bq  are  the  integration  constants. 

In  order  to  determine  U^,  the  general  solution  of  the  second 
of  the  differential  Equation  37  may  be  obtained  from 

dR2  R  “ 

2  2  2 
where  R  =  N  r 

The  general  solution  of  Equation  50  may  be  written  as 


(49) 


(50) 

(51) 


=  A1I1(R)  +  BiK1(R)  (52) 

The  boundary  conditions  in  terms  of  the  new  independent  variable  x  de¬ 
fined  by  Equation  46  are 


uo(Xi)  =  0 

Uq(x2)  =  U1(Ne) 

du  H/i  n  dU 
_ o  _  N(1  -  e)  _ i 

“  dx  2/3  ,2/3  dR 

al  b 

U^o)  =  0 


(53a) 

(53b) 

(53c) 

(53d) 
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where 


2/3  ,2/3 

X1  =  ai  b 


H 


2/3  ,2/3  ,.2 
x2  =  a^  b  /b 

From  boundary  conditions  53d,  we  have 

AiI1(o)  +  B^Co)  -  0 


(54) 


Therefore, 


B, 


0 


(55) 


Substituting  the  general  solutions  given  by  Equations  49  and  52  and  the 
condition  given  by  Equation  55  into  the  boundary  conditions  53a,  53b  and 
53c,  we  obtain 

A  AI(x- )  +  B  BI(x-)  =  0 

O  i  O  1 


A  AI(x-)  +  B  BI(X0)  =  A. I, (Ne) 
o  o  lx 


(56) 

(57) 


2/3b2/3  p 

-  - -  A  AI' 

al  [_° 


A.  I  (Ne) 
i  o 


(x,)  +  B  BI 
l  o 


I1(Ne)"*| 
"  ~ Ne 


’  (x2}j 


(58) 


where  the  primes  refer  to  differentation  with  respect  to  x.  A  ,  B  and 

‘  o  o 

A^  may  be  eliminated  from  Equations  56,  57  and  58  to  give  the  charac¬ 
teristic  equation  as 


-1/2 


AI'(x2)  BI(Xl)  -  BI'(x2)  AI(xx) 
^AKx^  BI(x2)  -  AI(x2)  BK^) 


I  (Ne) 
o  1 

Ix(Ne)  "  Ne 


From  the  definition  of  x^,  and  x2  in  Equation  54,  we  obtain 

1/2 


X1  =  X2  ~  "  e^x2 


(59) 


(60) 
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For  given  N  and  e,  the  values  of  x^  and  which  satisfy  Equations  59 

and  60  simultaneously  are  the  desired  values.  However,  not  all  pairs 

2 

of  x^  and  x^  lead  to  admissible  values  of  K  .  From  Equation  61  we  ob¬ 
serve  that  since 


K 


2e 


2e 


1  -  e 


X2  “  X1 


2  X2 


3/2 


(61) 


N(1  -  e)  ‘ 

2 

e^id  since  K  is  a  real  number,  first  of  all  must  be  a  positive  number. 
Secondly,  (x^  -  x^)  must  also  be  positive*.  Thus,  of  all  the  values  which 
satisfy  Equations  59  and  60,  we  select  only  those  pairs  of  x^  and  x^ 
values  which  satisfy  the  conditions  that  x^  >  0  and  that  x^  >  x^. 


Case  iit  The  Formal  Solution  for  m  =  0 


In  case  we  wish  to  consider  the  complete  range  of  e  in  1  >  e  >  0, 
we  must  consider  Equation  37  without  making  any  approximation.  We  then 
have 


where 


d2U  ,  dU 

_ o  _1  _ o 

,  2  r  dr 

dr 


v  .  2 

~2  +  “ 
r 


U  =  0 
o 


v2  =  1  + 


2  2 
AN  eZ 


2  -  XT2 

a  =  N  - 


K2(l  -  e2)2 
AN2 

2\2 


K  (l  -  e  y 

Equation  62  would  need  to  be  considered  with  the  second  of  Equation  31 
and  boundary  conditions  39.  The  general  solution  of  Equation  62  is 


(62) 


(63) 

(64) 


Uo  *  *v(ar)  (65) 

where  £  is  a  general  cylinder  function  of  order  v.  The  boundary  condi¬ 
tions  require  A  to  vanish  at  r  ■  1.  The  required  solution  may  be 
expressed  as 
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Uq  -  M  (J_v(a)  Jv(ar)  -  Jy(x)  J_v(ar)]  (66) 

where  M  is  a  constant. 

Defined  in  this  manner,  Uq  clearly  vanishes  at  r  =  1.  We  have  further 
at  r  =  e. 


M[J_v(a)  Jv(ac)  -  Jv(°0  J_v(ac)]  55  AiI1(Nei)  (67) 

Moreover, 

diJ 

o 


The  integration  constants  M  and  A^;  may  be  eliminated  from  Equations  67 

and  68  to  yield  the  characteristic  equation.  For  arbitrarily  assigned 

2 

v,  one  may  compute  K  from  the  characteristic  equation  and  Equations  63 
and  64  for  given  N  and  e. 


(Ne) 

VNe)  "  -5T- 


(68) 


THE  GENERAL  CASE 


Consistent  with  boundedness  of  the  solutions  of  and  R^  at 
the  origin  Equation  31  has  the  following  solution 

R,  =  AI  (Nr) 
l  m 

U  =  i  R'.  =  A  -  I  +  NI 
i  K  i  [r  m  mt-1) 

where  A  is  any  arbitrary  constant.  Equation  28  may  be  rewritten  as 


(69) 

(70) 


(71) 


(72) 
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The  boundary  conditions  become 

r-  — \ 

°o<«>  •  K 

f  I  (Ne)  +  NI  (Ne) 

e  m  m  +  l 

(73) 

R  (e)  =  AI  (Ne) 
o  m 

(74) 

U  (1)  »  0 

0 

, 

(75) 

where  Rq  is  redefined  as  -iR  of  equations  prior  to  Equation  69.  The 
eignevalue  problem  now  is  the  following:  find  those  values  of  K  for 
which  the  Equations  71  and  72,  and  the  boundary  conditions  73,  74,  and  75 
are  satisfied. 

Since  Equations  71  and  72  are  homogeneous,  any  constant  times 
any  solution  that  satisfies  the  boundary  conditions  will  also  be  a  so¬ 
lution.  Thus,  in  the  eigenvalue  problem,  we  can  arbitrarily  assign  any 
value  to  the  constant  A.  For  convenience,  we  set  it  equal  to  unity. 

NUMERICAL  SOLUTION 

Since  Equations  71  and  72  are  too  cumbersome  for  closed  solu¬ 
tions,  we  employ  numerical  integration  methods.  The  numerical  procedure 
for  solving  the  eigenvalue  problem  is  as  follows:  for  an  assumed  value 
of  K,  we  can  compute  the  values  of  UQ(e)  and  RQ(e)  for  a  given  set  of 
values  of  m,  e,  and  N.  Using  these  as  the  initial  values,  we  can  inte¬ 
grate  Equations  71  and  72  by  numerical  methods  to  yield  Uo(l).  If  this 
value  is  zero,  then  our  choice  of  K  is  indeed  the  right  one.  If  not, 
we  change  our  assumed  value  of  k  by  a  selected  increment  and  recompute. 

In  practice,  since  it  is  uneconomical,  computer  time-wise,  to  calculate 
a  precise  zero  value  for  Uq(1),  we  arbitrarily  say  that  if  the  absolute 
value  of  UQ(1)  is  less  than  a  prescribed  positive  small  number,  then  our 
choice  of  the  k  value  is  the  right  one. 

For  negative  values  of  m,  Equations  71  and  72  have  a  singularity 
depending  on  whether  or  not  the  quantity  1  +  (1  -  e“)  K /m  is  real  or 
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complex.  If  the  quantity  is  real,  the  singularity  is  at  a  raduis  given 
by 

'singular  '  »  '  fl  *  «  ‘  *2>  C/" 

At  this  radius,  the  derivates  of  Uq  and  Rq  cannot  be  computed  from  Equa¬ 
tions  71  and  72.  However,  from  physical  considerations, , we  require  that 

U  and  R  be  continuous  at  this  point, 
o  o 


When  a  singularity  exists,  we  employ  the  following  procedure: 

We  divide  the  two  zones,  e  <  r  <  r  .  and  r  .  <  r  <  1  into  N.  and  N0 

—  —  sing  sing  -  —  12 

intervals,  where  and  are  selected  to  be  such  that  the  increments 

in  radius  for  an  interval  is  approximately  equal  in  both  zones.  Starting 

from  r  ■  e,  we  integrate  the  Equations  71  and  72  till  we  reach  the  last 

interval  of  e  <  r  <  r  .  .  We  subdivide  this  into  four  smaller  inter- 

-  -  sing 

vals  and  integrate  the  equations  to  obtain  the  values  of  U^,  Rq  and  their 
derivatives  at  the  three  intermediate  points  shown  in  Figure  1.  Using 
these  values,  by  linear  extrapolation,  we  estimate  the  values  of  Uq,  Rq 
and  their  derivatives  at  the  singularity.  Since  we  cannot  use  Equations  71 
and  72  at  the  singularity,  using  the  estimated  derivatives,  we  integrate 
numerically  to  the  one-fourth  point  of  the  next  interval.  From  here  on, 
we  revert  to  the  conventional  integration  procedure  using  Equations  71 
and  72. 


The  numerical  integration  procedure  used  in  the  computer  pro¬ 
gram  for  calculating  the  eigenvalues  is  the  familiar  Runge-Kutta-Gill 
method. 

DESCRIPTION  OF  THE  COMPUTER  PROGRAM 


For  computing  an  eigenvalue,  the  input  data  to  the  program  is 
supplied  by  two  cards.  The  input  map  for  the  program,  describing  the 
breakdown  of  the  entries  in  the  cards  are  shown  in  Table  1.  The  third 
I-field  on  the  first  card  of  each  set  was  originally  intended  to  supply 
the  value  of  n.  Since,  this  was  subsequently  generalized  to  be  a 
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Table  1 
INPUT  MAP 


Variable 

Definition 

Card 

No.  1 

M 

Run  Number 

MM 

-m 

NN 

Enter  zeros 

NR 

Number  of  integration 
intervals  <_  500 

IGEN 

Number  of  Eigenvalues 
required 

NSTOP 

If  last  set  of  data  being 
read,  should  be  entered  non¬ 
zero.  Otherwise,  enter  zero 

* 

LC 

10  Ne  _<  50 

Card 

No.  2 

PLAO 

Initial  trial  value  of  Eigen 
value 

DPLA 

Increment  of  Eigenvalue 

TRUNC 

Permissible  truncation 
difference 

CRAD 

e 
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factor  in  N  (CAPN  in  the  program),  it  may  be  filled  in  with  a  zero  or 
any  other  number.  This  number  wil?  appear  on  the  printout  as  the  value 
of  n  (N  in  the  program).  Since  it  has  no  other  effect  on  the  program, 
it  may  be  ignored  for  all  other  purposes. 

The  computer  program,  written  in  Fortran  IV,  is  given  in  the 
Appendix:  the  first  subroutine,  NING(NRl)  provides  the  ^.unge-Kutta-Gill 
integration  method.  It  also  calls  a  subroutine,  DER(Y,  FM,  C0N,  CRAD, 

FK,  AK)  which  provides  the  derivatives  of  Uq  and  Rq  as  given  by  Equa¬ 
tions  71  and  72.  Subroutine  NUESTI,  called  by  the  main  program,  provides 
a  method  for  incrementing  the  eigenvalue  on  the  basis  of  current  deriva¬ 
tion  in  Uq(1)  from  the  specified  truncation  error.  Subroutine  RKGC0N 
(AIN,  BIN,  CIN)  supplies  the  RKG  constants  to  the  main  program.  Sub¬ 
routine  BESC0N(AI,  BI) ,  called  by  the  main  program  supplies  the  values 
of  the  modified  Bessel  functions.  If  the  program  is  to  be  run  for  any 
values  of  m  other  than  -1,  this  subroutine  will  need  to  be  changed  to 
supply  the  modified  Bessel  functions  of  the  appropriate  order. 
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RESULTS 

Table  2  shows  the  effect  of  various  step  sizes  on  the  computed 
eigenvalues.  Truncation  is  of  10  **,  initial  eigenvalue  is  0.01,  eigen¬ 
value  increment  is  .05. 

Table  2 

EFFECT  OF  STEP  SIZE  ON  COMPUTED  EIGENVALUES 


Lowest  eigenvalue  K 
Ne  =  .2 

Total  No.  of  b/a  =  e  =  0.45 

Steps  _  m  =  -1 


100 

.76079984 

200 

.75407226 

300 

.75447285 

400 

.75445469 

500 

.75479767 

for  Lowest  eigenvalue  K  for 
Ne  =  .2 
b/a  =  e  =  .1 
m  =  -1 


.35125716 

.35046438 

.35032046 

.35025621 

.35021697 


The  results  show  that  if  200  or  more  total  number  of  steps  are  employed, 
the  variation  is  in  the  fourth  significant  digit. 

Table  3  lists  the  results  for  various  Ne  and  e  values  computed 
with  200  integration  steps,  an  initial  eigenvalue  of  0.01,  and  an  eigen¬ 
value  increment  of  0.05.  With  this  sweep  procedure,  the  listed  eigen¬ 
values  were  found  to  be  the  lowest.  It  is  possible  that  if  the  eigenvalue 
increment  is  lowered  from  0.05  to  some  smaller  value,  other,  ever  lower, 
eigenvalues  may  be  discovered;  but  this  is  considered  to  be  extremely  un¬ 


likely. 


In  Stewartson’s  notation,  we  have 

C  7T 

a(2j  +  1)  “  2N 
b/a  =  e 

_ t  =  K _ 

Stewartson,  K. ,  On  The  Stability  of  a  Spinning  Top  Containing 
Liquid ,  J.  Fluid  Mech.,  5,  1959,  pp.  577-592. 
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APPENDIX 


A 


REOU I 

main 


04/01/66 


EFn  source  statement 


I fN ( S  > 


MAIN  PROGRAM  FOR  ThE  DETERMINATION  OF  EIGENVALUES  31G-B2294-01 

nCOMMQN  FM,FN,CON,DEL.NR.X(1,501>.Z(1#501).FRE8U(3)#DET(3),AIN(4), 
1BIN(4) ,CIN(4)  .CRAD.DPl*. I DET , I LU , KLUE , UCO , RCO , R AD  I  US ( 50 1 ) 

DIMENSION  A  I  <50  )  # B I <  50  > 

NSTOP»0 

CALL  RKGC0N(AIN,8IN«CIN) 

CALL  BESCON(AI.BI) 

20  IF(NSTOP.NE.O)STOP 

READ<5, 10 )M .mm.nn.nr, igen,nstop.lc,plao.dpla,trunc,crad 
10  FORMAT(7I6/3E12.0.F12.6) 

MMs-MM 
FM  =  MM 
FNsNN 
FNRsNR 

FREQU ( 1 ) sPLAO 
BNC=LC 

fnc=bnc/io • 

CAPN*FNC/CRAD 

C0N=CAPN»*2 

a!mi«ai (lc)*exp(Fnc> 

a!M2=BI <LC)*EXP(FNC> 

IOET»l 
K«K  =  1 

WR I Tt  t  6  »  85  ) 

fi5  FORMAT { 1H1  35HCURRENT  VALUES  Of  INTEGRATION  DATA  //) 
wRITt(6.66)CAPN,LC 

86  FORMAT ( 3X7HCAPN  «  Fi0.4,5X5hLC  ■  12) 

DO  30  I  G«1 » I  GEN 
I  COUNT *1 

70  UC0= (FM*AIM1/CRAD*CAPN*AIM2)/FREQU(1) 
rco=aimi 

RADIUS(1)=CRAD 
X  < 1 > 1 1 *uco 
Z<1»1)*RC0 
CALL  NJNG(NRl) 

DET(1)«X(1,NR*1) 

WRITE<6«90)M,IG»FREOU(1>.X(1iNR*1) 

90  FORMAT (  //  10X  13HRUN  NUM 

tBER  a  13.  5X  5HIG  »  12,  5X  16HTRIAL  E I GENV .  «  E15.6,  5X  9H  U  < 1 > 
2=  E32.0) 

I F ( IDET-2)170. 170.160 
170  Dl  *  DET  < 1 ) 

AL  =  FREQU(l) 

GO  TO  190 

180  If(AL»FREQU(l) J170.170.190 
190  ABET  *  ABS  ( DET ( 1 )  ) 

IP(ADET»TRUNC)80»80»100 
100  CALL  NUESTI 

GO  TO  (80,110,  20),KLUE 
110  I OET  »  I  DE T *  1 

I  COUNT  *  I  COUNT ♦ 1 
IF( ICOUNT-50)70,70,  20 
80  WR  I  Tfc ( 6  » 1 ) 

1  FORMAT ( 1H1  41X  32  I  EIGENVALUES  and  EIGENFUNCTIONS  ) 


G0000013 


00000014 

00000016 


00000016 

00000017 


00000023 

00000024 

00000025 

00000026 

00000027 

00000028 

00000029 

00000030 

00000031 

00000032 

00000035 

00000036 

00000037 


00000040 


A1 


04/01/66 


ME  DU  I 
main 


EFN  SOURCE  STATEMENT  -  J  FN  <  S )  - 


Wl’!Tb(6,2)M,HH,NN,  IG,FREQU(1)  .CAPN.CRAD 
2  FORMAT! //4dx  10HR  IN  NC.  =  13///  7X  6H  M  *  1 3 .  5X  8h  N  *  1 3  > 

1  6X  12,  2X  13HP1IGENVALU6  *  E15.8.9X6HCAPN  *  Fi2.6,5x4HC  *  F6.2) 
WM  j  T  t  <  6 , 3  ) 

:  FORMAT ( ///19X  2H  R  4X  1AH  u  MODE  SHAPE  4X  14H  R  MODE  SHAPE  18X  2H 
3R  4X  1 4 h  u  MODE  S  Ur'E  3X  15H  R  MODE  SHAPE  //> 


Kp  A*  1 
N«R=NR/2 
qM  40  1*1, NRR 
IF  (i=6t*KRA)  60,60,50 
60  KpA  1  Kr>«  +  1 
Wp I Tt ( 6 , X  ) 


00000047 

00000048 

00000049 

00000051 


WRITc(6,2)M,MM,NN. I3,fR6QU(1),CAPN,CRAD 
WR I Tfc ( 6 , 3  ) 


50  K»NRk*I  00000055 

WR I Tt ( 6 , 130 )  RADIUS!  1  ),X(1, I ),Z<1. I > , RAD  I  US < K > , X < 1 , K ) , Z < 1 4 < > 

130  FORMAT ( F20 • 4 ,  3X  E15.8,  3X  El5.8,F2Q.4,  3X  E15.8,  3X  E15.8)  00000057 

IF(I- NR R)  40,140,40  00000058 

140  K  =  **1  00000059 

WRITe<6,150>  RAD  I  US ( K ) , X ( 1 , K  ) , 2 ( 1 ,  K > 

153  F')RmaT(40x  F20 . 4 , 3  x  El5,8,  3X  E15.8) 

40  CONTlNUt  00000062 

FREOUC2)=AL  00000063 

DFT(2>  =DL  00000064 

FREQO(l)=AL+DPLA  00000065 

30  IDETS2  00000066 

G'1'  TO  20  00000067 

END  000 00068 
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REDD! 

GOODY  -  EFN  SOURCE  STATEMENT  -  JFN(S)  - 


I 

I 

r 


1  2 

! 


3 


1 

i 


li 

I 

I 

I 


50 

100 


1000 

} 

1 


i 


r 


SUBROUTINE  NINQ(NRl) 

OCOMMON  FM|FN,CON|DEL»NR»X(1#501)iT(1i501)«FREOU<3)#DET<3)*A!N(4) 
1B!N(4) .CIN(4),CRA0,DPI A. !DET, ILU,YLUE,UC0,RC0.RADIUS(501) 
dimension  ak<3>,y<3>,o<3>,am(3>,bm<3> 

FK«FrEQU(1) 

Y ( 1 ) *CRAD 
Y ( 2 ) *UCO 
Y ( 3 ) s  RCO 
Q(1>=0, 

Q<  2 ) *0  . 

Q( 3  >  *0 • 

AK(1)«1, 

B!NG«1.*(1,-(CRAD**2)  )*FK/FM 
IF  ( 9 1 NG ) 3 1 3 « 1 
RSINQ*CRAD/SORT(BINQ) 

WR[TE(6<2)  RSING 

FORMAT ( //20X  14HR(SINGULAR)  *  F9.fl 

GO  TO  4 

WRITEI6.5) 

FORMAT(10X3#HSINGULARITY  OUTSIDE  INTEGRATION  DOMAIN  1 

BNR»NR 

NR1«8NR*<RS!NG-CRAD)/<1.-CRAD> 

NRMsNR-50 

IF(NR1.LT'.50)NR1*50 

if(nri.gt;nrm)nri*nrm 

NR2*NR-NR1 

FNRlsNRl 

FNR2=NR2 

DELa(RSlNG-CRAD>/fNRl 

DELl»<l,-RS!Nfl)/FNR2 

DFLTA»DEL/4. 

nAB«NR1*2 

DO  1000  UpI.NAB 

IF(L ,GE.NR1»DEL*DELTA 

DO  100  JJ«li 4 

CALL  DER( Y.FM.CON.CRAD.FK, AX) 

DO  50  1*1,3 

AIKN  sA!N(JJ)*<AK(!)-BIN( JJ)*0(  Ill 

YU)  *  Y(  I  )*DEL*AIKN 

om  *  0  <  I )  *3  i  *A  I  KN  -CIN(  JJ)«AK(  I  » 

CONTINUE 

RAD!  JS(L*1)«YU) 

X ( 1 , L*1 >  * Y<  2 ) 

Z  < 1 • L*l 1 • Y ( 3 1 

CALL  DER( Y»FM,C0N.CRAD.FK. ak  ) 

B  U  *  A  K  (  2 1 
BR*AK(3> 

Y(l)aRADIUS(NAB) 

CALL  DER( Y,FM,CON,CRAD»FK, AK) 

A  LJ  S  A  <  (  2  > 

AR*AK(3) 

AK(2)*2,*GU-AU 
AK ( 3 ) *2 . *IR»AR 
Y<1)»RS!NG 

Y(2I»2.*X(1,NA8*1)-X(1,NA§) 


05/02/66 


A3 


05/02/66 


REPDI 

GOODY  -  EFN  SOURCE  STATEMENT 

Y(3)«2.*ZU»NA8*1)-Z<1,NAB) 

RAD! JS<NR1*1)»Y(1) 

X<1,NR1*1>«Y<2) 

Z!1.M»1*1>«Y(S) 

DEL'DELl/4, 

DO  80  L ■ 1 1 4 

DO  80  JJ* 1  •  4 

!F(L.EQ,1>G0  TO  81 

CALL  DER( Y,FM,CON,CRAn»FK,AK) 

«1  DO  82  I  * 1 »3 

A ! KN  «AlN| JJ)*(AK( ! > -B I N < J J ) *0 < I )) 

Y( T )  •  Y< | )*DEL«AIKN 
82  0  (  i >  *  0 ( ?  >  *3 . *A I KN  -CIN( JJ)*AK( I > 

80  CONTINUE 

RADI JS(NR1+2)»Y(1) 

Xtl*NRl*2)*Y(2) 

Zfl.NRl*2)«Y(3) 

DEL»DEL1 

DO  900  L*NAR , NR 

IF(L,E0.NR)DELS1«"Y<1) 

DO  9o  JJ»1»4 

CALI  DER( Y,Fm,CON,CRAD,Fk, AK ) 

DO  60  I >1 1  3 

A  I  KN  =  A  I  N  *  JJ>*(AK(  I » -s-B  I N  (  JJ)*0(  I  >> 

Yin  *  y<i)*del*aikn 
60  OU)  »  0  (  I  )  *3  ,  **  I  KN  -CJN(  JJ)*AK(  n 
90  CONTINUE 

RAD  I  US ( L*1 ) *  Y 1 1 > 

X(1»l*1)*Y(2) 

900  Z(1.l*1)»Y(3) 
return 
END 


IFN(S)  - 


09/F9/66 


REDD! 

DERIVE  -  EFN  SOURCE  STATEMENT  -  IFN<$)  • 


SUBROUTINE  DER( Y#FM# CON, CR AD i FK ,  AK) 

DIMENSION  Y(3)»AK(3> 

R?s1./(Y(1)«Y(U) 

ALPHA=FM#FM»R2*CQN 

BETAs2./(l,-CRAn*CRAD) 

OMEGAs .5#(lB-CRA0*CRAn#R2 )#BET A 
SIGMAsFK4FM*OMEGA 

AK(2)=<ALPHA*Y(3)*FM#BETA*Y<2)/Y(i  ) ) /S I GM A- Y ( 2 ) / Y ( 1 ) 

AK<3)  =  -2,#0MEGA*(BETA*Y(2)*FM*Y(3>/YU)  > /S I GMA*S I GMA«y < 2) 

RFTURN 

END 


1 


I 


1 


1 

l< 

1 


1  J 

v 


41 
31 
31 
3  1 
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KFrmi 

PPI  AR  -  EFN  SOURCE  STATEMENT  -  I FN C  S I 


SUBROUTINE  NIJFSTT 

OCCMMIIN  F  M, F  N. CUN. OFl  .NR.XI  1 . 50 1 ) • 7  I  1 , 50 1 ) , FREQU ( 3 ) * DET ( 3 ) , A  IN ( 4  I , 
1RINI41,CIN(  4 )  .CRAD.DPl  A,  10FT,  ILU,KLUF,UCO  .RCO.RAOIUSI  50  1) 
n:i)IIRi  F  PRECISION  F3?l,F213,FD,AFl,F23,F3l,Fl2,OiFi,D2F2,D3F3*AP2, 


I  c  I  42 • AF3.QD. CC, SR  1 ,FF 1 , FF2 
I F  (  l OF  T- 2 )  10,  20  •  30 

10  FRFOUI  ?>=FRFQIJ(  1  ) 

F R F 0 U  (  1  )  =FRFOtJ(  2  )  4-DPL  A 
CFTI 2)=nFT( 1 ) 

NCI  UF=  3 
GO  TO  40 

20  CDF  I  =  OF  T (  1 )*  OFTI 21 
SCIDF  =  3 

I F ( D  0  F  T I  SO  «  SO,  60 

60  II  U  *  7 
GU  TO  70 
SO  III)  *  1 

7P  FRFOUI 3)=FRF0UI 2) 

C F T ( 3 )  =  OFT ( 2  $ 

FRFOUI 21  =FRFGU( 1 ) 

OFTI  2)  =  OFT  I  1 ) 

OS  Gfi  TO  I  90 •  1  on )  •  II  u 

OO  FRFOUI  1  ) =. 5* I  IFRFOUI  3  )  4-FPFQUC  2  )  I  — (  (FR  OUI  3 1 -FRF GUI  2)  I/IOETI  31- 
x  OFTI  2m*(DFTI3U0FTI  211  ) 

GO  TO  40 

100  FRFOUI  1  ) -FRFOUI 21+0 PL  A 

r;n  to  40 

'  i.n  TOC  1  10,1  201,  ILU 
1  10  OOF T  =  OFTI 1  >*OFT ( 2  I 
IhlDDFTI 130, 130,140 

1  30  NCI  UF  =1 
GU  TO  1  SO 

140  nCFTaflFTI  1)*DFTni 
Nfl  UF  =  2 
I  FI  OOFT 1150,1  SO, 70 
1  SO  f FI  FRFOUI 31  I  165,  165.  300 
165  rn  TDI 17S,1 RSl.NCLUF 
1  7S  FRFOUI 3  I =FR  FQU(  1 ) 

OFT  I  3 ) =OFT I  I  1 

go  rn  95 

1  PS  FRFOUI  2I=FRFQIJ(  1  I 

nm  ?  i  =dft 1 1 1 

GP  TO  os 

300  0  F  1  2  =  A  R  S  I  IFPFOIJI  ll-FRFOlM  2  I) /FRFQUI1)  1 
on  V-ARS  I  IFREOHI  1  l-FRE  OUC  3  )  1  /FREOlJ(l) ) 

OF 2 3  =  A RS  I IFRFOUI21-FREOUI35 )/FRE0UI211 
IF! OF  1 2-4  IF- 04)312, 312. 313 
111  1  FI  OF l 3-  . 1F-04)312, 312.314 
314  IFI0F23-.1F-041T12.312.3U 
31?  <;o  rru  16S,soi,  ilu 
311  01  *  OFT  I  1  I 
C2  ^  OFTI 21 
03  =  OFTI 3) 

FI  =  f  PFOUI  1  I 
F  ?  =  FRFOUI 2) 


00000i?3 

0O000124 


00000126 
00000127 
0OOO0128 
00000129 
00000130 
pnnOO 131 
00000132 
00000133 
00000134 
0O0001 3S 
00000136 
O0000137 
00000138 
noonoi 39 


O00C0141 

00000143 
00000144 
O000014S 
00C00146 
0OC00147 
O000014B 
00000149 
00O001S0 
00000151 
00000 152 
00000153 
00000154 
00O00155 
00000156 
00000157 
00000158 
00000159 
00000160 
00000161 
00^00162 
00000163 
00000164 
00000165 
00000166 
00000167 
00000168 
00000169 
00000170 
00000171 
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rfooi 

Pfll  AR 


-  FFN  SOURCE  STATEMENT  -  IFN(S)  - 


FT  *  FRFGUC3I 

F37T  *TF3-T71/F1 

FIT?  =(F1-F3)/F7 

F ?l 3  *(F?~Ftl/F3 

FD  *  F321  ♦  FI  3?  ♦  F ? 1 3 

AF1  *ini#F3?l  ♦  02*F 13?  ♦  D3*F713I/FD 

F73  =  F7/F3 

FT!  =  F3/FT 

FI?  =  Fl/F? 

01F1  =  Ol/Fl 
n?F?  =  n?/F? 

D3F3  *  D3/F3 

AF?  *  U)iFl#<F?3-1.0/F23l*D2F2*(F3l-1.0/F3i)*03F3*(F12-l 
1  FD 

AF3  *  <D1F1*(  1.0/F2-1.0/F3I*D2F2*U„0/F3-1.0/F1I*D3F3*'<  1 
1/F21l/Fn 
OP  =  AF1/AF3 
r.C  =  0.5*AF?/AF3 

r.cn«rx**?-o.) 

IFTncm317.312.397 
SRI  =  SORT  ( CCD ! 

FF1  =  -CC+SR1 
FF2  =  -CO  -SP1 
GO  TOC 160.40CI# It  U 
*  1UJ  »  1 
FRF0TTT3T  =  FRFC 
DF  (  il  =  DFT i ^  I 
FRFGIJ(7)  =  FRFOlim 
OFT  C  7 1  =  DFT( II 


T7TT 


IFC  FF 1-FRFOUI 31  1500*  500*600 
FRFOUT'T 1 =FF7 
GO  TO  40 

IFf  FF1— FRFOIJC  71 1700*  700*  500 
FRFQUf 1  I  =FF1 
GO  TO  40 

GO  T0(  1  70.  I  80  I • NCLUF 
FRFOirt  if=FRFOTrT  IT 
OFT  C  3  I  =0FT (II 
GO  TO  BO 

FRFOIM  ?)=FRFQU( 1 1 
DFT C  7 1  =OETf 1 1 
GO  TO  BO 

TTTTFT  =  OfTM  1*ITFT1 77 
IFIOOFTMSO# ISO#  200 
NCLUF  =  7 

IFCFRFQUI3I ISO#  50*  300 
NCI  UF  =  3 
GO  TO  70 

xo  mrTiTT#  77rr#7wr#  ncxttf 

OFR  A  =  ABS  C  FRFOUI 1 l-FRFOUC  311 
GO  TO  740 

OFR  A  =  ARS  (FRFOUI l l-FREOUl  211 
TFC DFRA-c 1 F— 071 250#  250#  230 
KiUF  =  l 

GTTTnT77STT#7TTT#  7  ROT  .TCIIJF 


02/22/66 


0°000172 
00000173 
00000174 
00000175 
00000176 
00000177 
00000178 
0000017? 
00000180 
00000181 
00000182 
00000183 
,0/F 12  1 1 /000001 84 
000001 B5 
,0/Fl-l. 000000186 
00000187 
00000188 
00C001 89 
00000190 
00000191 
00000192 
00000193 
"0000194 
00000195 
00000196 
00000197 
00000198 
00000199 
00000200 
0 0000201 
00000202 
00000203 
00000204 
00000205 
00000206 
00000207 
00000208 
00000209 
00000210 
00000211 
00000212 
00000213 
00000214 
00000215 
00000216 
00000217 
00000218 
00000219 
00000220 
00000221 
00000222 
00000223 
00000224 
00000225 
00000226 
0"000227 
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PDI  AR  -  FFN  SnURCF  STATEMENT  -  IFN<S»  - 

260  FRFOUI  1  I  =  FRFOlim  00000228 

On  TO  280  00000229 

'  1  FRFOlim  =  FRFQUm  00000230 

or  TO  PRO  00000231 

230  K 1  II F  =  2  00000232 

280  CflNTINUF  00000233 

RFTIJRN  00000234 

FN(1  00000235 
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RUNqe  -  FFN  SOURCE  STATEMENT  -  irNfS)  - 


subroutine  rkrconi ain»bin»cin> 

DIMENSION  A  I N  C  4  > »RIN<4) »CIN<4> 
AlN(l)sl./2. 

SRT=SORT( AIN(l)  ) 

A  I N ( 2 )  =  1 . -SRT 
AjN(3)sl.+SRT 
AlN<4)sl,/6. 

B l N  < 1 ) =2 . 

BIN(2)=1 . 

8  I N ( 3  )  si . 

B l N ( 4 ) s2 • 

C I N ( 1 )sAlN(l  ) 

CIN(2)sAIN(2) 

CIN(3)=AIN<3) 

CIN(«)=A1N<1 ) 

RF  T  URN 
END 
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SUBROUTINE  HESCON(  A  I  ,  R  J  ) 
DIMENSION  AI(50) .RJ(50) 

A  f  < 1 )  *0,0452904 
AI  (?)  *0.0822831 
M(3)  *0.1123775 
A  T ( 4  )  sO,  136763? 

A  j  <  5  J  *0  *  156420® 

A  I ( 6 )  *0.1721644 
AI (7)  *0,1046699 
A  I  ( ft )  *0.1  944907 
A  I  <  9  )  =0.2021165 
AT (  1  0  1  =  0 . 2079104 
A! ( 11 >  *0 , 2122016 
AI  (12»=0. 2152568 
A  1 ( 13 ) *0 . 2172976 
AI (141*0.2185076 
AT  <15>=0. 2190694 
AI (16)*0, 2190195 
AI (1/1=0.2185528 
AT (10)=O. 2177263 
AI ( 19 >=0 . 2166120 
AI (?0>=0. 2152693 
AT (21>=0, 2137478 
AI (221=0.2120877 
AT (231=0,2103230 
AT (241=0.2084811' 

A! (251=0.2065846 
A  K  26>  =  0 . 2046523 
AI (271=0, 2026990 
AT (201=0,2007374 
AT(291=0. 1987773 
A  I  (  30 1 =0 . 1960267 
A[(31>=0, 1948921 
A  I  (  321  =  0  .  1  929706 
AI  (331=0,1  910902 
AT(34)=0, 1892299 
A  I ( 35 ) *0 , 1  873999 
AI  (36  1=0,1856022 
AI(37)=0, 1838379 
AI  (381=0,1821076 
AI  (391=0.1004119 
AI  (4O'=O,17fl7508 
AI  (41  1  =0  ,  1  771244 
AI  (421=0,1  755325 
AI  (431*0.1739746 
AI  (441=0.1724502 
AI  (451=0.1709588 
A  I  (46  )  =0 . 1694997 
AI  (47  1*0.1680723 
AT (481=0,1666757 
AI <491=0.1653093 
AI <501=0,1639723 
BKl)  =,90  71009 
B I ( 2 )  *.8269385 
B I ( 3 )  =.7575806 


iSljftM 

actually 

ISUB1 


I  S(JBM*1 

actually 
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SOURCE  STATEMENT 


IFN(S) 


02/09/66 


RtIUlI 

GOOF  4 


EFN 


8114)  a,6974022  ISUBZERO 

B 1 ( 5 )  a , 6450 

8 1 ( 6  )  s,5993272 

81(7)  8.5593055 

81(8)  8,5241489 

B 1  ( 9  )  b,493163o 

B1  <10>s. 4657596 

8Kll)a, 4414404 

81  (12>  =  , 4197821 

81  (I3)s, 4004249 

81(14)8,3830625 

B1  (15)  =  , 3674336 

81  (16)8.3533150 

81(17)8,3405157 

81(18)8,3288719 

B1 (19)8,3182432 

81 (20)8.3085083 

01(21)8,2995631 

8K22)  =  .29l3l73 

B1 (23)s,2836930 

81(24)8,2766223 

Bl(25)s,2700464 

81(26)8,2639140 

B 1  ( ^ /  )  s , 2581801 

B!  (28)8.2528055 

81  (29)8,2477557 

B!(30)s.2430003' 

BJ  (3l)s, 2385126 

01(32)8.2342688 

81(33)8,2302480 

B|(34)s,2264314 

81(35)8,2226024 

HI  (36)8.2193462 

61(37)5.2160494 

81 (38)8.2129001 

81 (39)8,2098675 

81 (40)8.2070019 

81(41)8,2042345 

81(42)8.2015774 

BI<43>8. 1990232 

01(44)8.1965656 

01 !45)s,19419e3 

81 (46)*.1919160 

81 (47)8,1697134 

81(46)8,1675662 

81 (49)8,1655300 

81(50)8.1035408 

HF  TURN 

END 
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